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Abstract 

We consider the Helmholtz problem in the context of the evolution of uniform initial distribution 
of a physical attribute in general porous media subject to a partially absorbing boundary condition. 
Its spectral property as a reflection of the boundary geometry has been widely exploited, such as 
in biological and geophysical applications. We consider the situation where the critical assump- 
tions which enable such applications break down. Specifically what are the consequences of an 
inhomogeneous absorption strength? By means of perturbation theory exact theoretical results, 
and numerical simulations on random sphere packs, we identify the regions of parameter space in 
which such inhomogeneity is important and those in which it is not. Our findings shed light on 
the issue that limits the mapping between the diffusion/relaxation spectrum and the underlying 
boundary geometry. 

PACS numbers: 89.90.+n,76.60.-k, 81.05.Rm,91.60.-x,82.56.,87 



The spatio-temporal evolution of the density of an attribute carried by diffusing agents, 
subject to a partially absorbing boundary, occurs in a variety of scientific context that ranges 
from NMR relaxometry in porous media|IJ[2] and bio- medicine [3], waves in a membranejl], 
to migration of genes and cultural influences [5]. In its minimal form, the problem is formu- 
lated as the Helmholtz problem[6], but in real systems, variations in the absorption strengths 
and local geometry complicate the matter. As a result, changes in part of its spectrum or 
phantom length scales, may appear. Empirical data from such systems, usually lacking the 
reference system with which to compare, beg the question: Is such a feature inherent in 
the nature of the underlying dynamics (with a clean boundary) or due to the haphazard 
elements on it? We are directly motivated by NMR relaxometry in porous media in which 
issues of such nature have been long standing. In this work, we consider the situation where 
the draining strength, p(r), varies from point to point on the boundary surface and probe 
how it intertwines with the boundary geometry. Although some studies exist on aspects of 
variable p(r)[3 El El EH], their systematic investigation is still lacking. In this letter, we de- 
velop a theoretical framework which incorporates both effects on an equal footing by treating 
the spatial fluctuations of p(r), namely 5p(r), as a perturbative parameter and identify key 
spectral features observable through numerics or experiments. The method is then applied 
to a simple problem, for which we obtain an exact solution for comparison. For realistic 
porous media and 5p(r) variations, we perform numerical simulations to determine bounds 
for the observable consequences. The result shows that the effect depends sensitively on 
the symmetry properties of both Sp(r) and the eigenmodes of the boundary- value problem 
associated with the uniform p. 

We start with a generic and widely studied problem: the evolution of a local density 
\l/(r, £) of a scalar attribute carried by entities diffusing (with diffusivity D) inside a general 
pore space (V p ) defined by the pore-matrix interface (£), which drains the attribute with 
a strength parametrized by the parameter po(> 0) upon contact with an entity. Defining 
the density operator as J = — _D(r)V, and TC = V • J, the local continuity of \I> leads to 
the classic Helmholtz equation, and in seeking its solution in terms of the eigenmodes (i.e. 
{<f)p(r)e~ x p t }), one arrives at : 

w#00 = K<t>» (!) 

where each mode is subject to the Robin's condition h ■ J0 p (r) = p o 0p(r) on the boundary. 
We superscript A° and 0° to indicate that it is for the case with a uniform p . From the 
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properties of H and po, it follows that X p are real, and > 0; orthonormal eigenmodes (j) p 
can be represented as real functions. The connection between the spectrum of 7i and the 
boundary shape (pore geometry) had been noted by Kac and others [TTI IT2"| [T5] and some of 
its consequences were widely exploited in a variety scientific disciplines [El [15J, [16l ITT] . In 
particular, the relaxation of polarized proton spins carried by fluid molecules in porous media 
such as rocks and biological samples is a pertinent example, as the relaxation is enhanced 
by the fluctuating dipole field near the interface. [HI [19] The evolution of the total attribute 
(magnetization in the case of NMR) M. = j v dr^(r,t), if we assume it to be initially 
uniform, then follows Ai(t) = ^ p e~ Xpt |a°| 2 where a° p = j y <pp(r)dr. In the case of a simple, 
closed boundary with a single defining length scale (such as the radius a in a spherical pore), 
this enhanced relaxation may be limited either by the strength of relaxation at the boundary 
or by the diffusivity D and a control parameter k = poa/D emerges to separate the regimes 
which have distinct spectral properties (i.e. the weights a p and values of the slow modes 
X p ,p — 0, 1, . . .). In the limit where k <C 1, it was observed that a ~ 1 for the slowest decay 
mode (p = 0), and also that A° = p S/V p (= po3/a for a sphere) directly proportional to 
the surface-to-volume ratio of the pore. In the other limit, faster modes generally gain in 
weight, and A[j ~ Aqo = Dir 2 /a 2 . In many situations, the relationship Aq ~ p/a is exploited 
to map the observed spectral distribution a p ,p — 0, 1, ... to a distribution of pore sizes. For 
this mapping to work, however, the pores should be isolated from each other (i.e. diffusive 
coupling [20] is negligible), each satisfying the condition k C 1, and finally p should be 
uniform. In an extended pore, the local variations in its geometry and the strength of 
the diffusive coupling make it no longer feasible to characterize the dynamics with a single 
length scale parameter a, and therefore the first two assumptions break down. Furthermore, 
in most real systems, p assumes a spatial variation (we will call it texture from now on). In 
addition to the fundamental issues as posed by Kac and others [HI [12], it is straightforward 
to show that, for a given distribution a p arising from a collection of pores distributed in sizes 
all satisfying k C 1, one can construct a collection of identical pores with an appropriately 
chosen range of p strengths assigned to each. 

In real life porous media, the pore space forms an extended, multiply connected manifold, 
with possible local variations in the connectivity of its constituents. It becomes necessary, 
then to employ notions more than pore sizes and throats for such systems, analogous to 
the way one progresses from atomic orbitals to the band theory and further on for disorder 
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in solid state physics. In this approach, a few properties of the eigenspectrum become key 
elements [2"T| 122] . In the following, we seek to quantify how inhomogeneous p(r) affects the 
eigen-spectrum and its experimental manifestation. 

Let us first generalize the Helmholtz problem Eq{l]by substituting p(r) = po + 8p(r) for 
po with the requirement §^ 5p(r)da = 0. Under this boundary condition, we pursue the new 
set of eigenmodes with eigenvalues X p . Using the self-adjointedness of Ti and Green's 
theorem, we obtain the following relationship: [2Tj 123] 

\ p = I p{v)\ < P P {v)\ 2 da + f D{r)\V<f> p (r)\ 2 dr (2) 
is Jv p 

which expresses all eigenvalues as a sum of surface integral and the diffusion-controlled 
volume integral, analogous to the energy of a particle in a potential well given in terms of 
the potential and the kinetic part. The character of each eigenmode can be understood in 
terms of the competition between these two components. A general observation can be made 
for fast modes that the second component dominates and X p for p > becomes progressively 
insensitive to p and 8 p. We define k p = by introducing a length scale parameter for each 
mode £ p = j^^^p ■ Applied to the p = mode, the criterion k <C or 3> 1 generalizes 
the earlier observations made for simple closed pore geometry[19] and is reminiscent of the 
A parameter for the electrical conductivity of pore filling fluid [21]. The spectral weight for 
the slowest mode ao is significantly weakened for Kq ^> 1, which leads to A4(t) with a multi- 
exponential characteristics that had invited the potentially misleading interpretation based 
on isolated pore size distributions. Instead, we derive a relationship that shows that this 
weight is directly related to the spatial fluctuation of the slowest eigenmode: 

\a \ 2 = l-V p ([ rfr|0 o (r)| 2 -| / rfr</» (r)| 2 ) . (3) 
v Jv p Jv p J 

Note that these rigorous relationships (Eq)2j [3]) apply to general boundary shape, and both 

uniform and inhomogeneous p. It is also straightforward to prove that slope of log.M(i) at 

early times should remain robust against the fluctuatons 8p(r), — lim t ^ ^l°g.M(t) — > poy 

but the range over which this is valid could be severely limited depending on the strength 

of \8p\. Many authors had considered the so-called mean lifetime r = ^ p a p /X P [USE], f° r 

which we obtain r = r — y- ^ dcrMo(r)5p(r)ito(r) where uq = lim s=0 J °° ^(r, t)e~ ts dt with 

\l/(r, i) being the local density under p(r), and similarly with m° and \l/o under the uniform 

Po- 
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The eigenvalue Ao and the spectral weight |ao| 2 of the slowest mode are the most accessible 
indicator for the change in the boundary condition p ~~ *• Po + ^p( r )- in the following, we 
therefore focus on the fractional shift in the slowest eigenmode = A °~ A ° which determines 
the long-time slope of logM.(t) vs. t. Figure [l] summarises schematically these general 
observations. We first derive a perturbative solution for ^jf for an arbitrary pore geometry 
and p(r) texture, and compare the result with an exact solution. Assuming that the complete 
eigenmodes {4> p } with A°'s are worked out already for the uniform p , we put the eigenmodes 
for the new boundary condition with p(r) as 

p (r) = c p (tj> p (r) + 5> re 0;(r)) + Q p (r) (4) 

q+p 

where c p is the normalization constant. We introduce the auxiliary function Q p (r) = 
p (r) — f dr'V(r,r')<p p (r') defined via the projection operator V(r, r') = J2 P 0p( r ')0p( r ) 
onto the Hilbert space spanned by the eigenmodes {<p p } of the uniform case. Note that 
formal inclusion of Q p (r) is necessary at this point to satisfy the new boundary condition 
as (f) p , if it were to be spanned by {4> p } alone, would satisfy the uniform p condition. Using 
the orthonormality of the complete sets {<f) p } and {4> p } respectively, it is straightforward to 
obtain a recursive equation for a pq : 



l a pr$Pqr + ^Pqpjy + $pq (5) 



> 9 r P 



where we define overlap integrals 5p qr y- = <f>^ da <p q 5 pep® and 8p qp y- = § d<j?,(f) q 5pQ p . Via 
iterative substitutions, we obtain the desired result in a systematic power expansion in 5p: 

k = K + y - £ + s to) + ( fi ) 

IT^P ' P P 

Defining / p (r) = — 4> q (r) §^ da(p q 6p(j) p , a representation of 6p(r) projected onto the 
Hilbert space spanned by {<p p }, we find that Q p should satisfy the inhomogeneous equation: 
(H - Ap)Qp(r) = f p (r) subject to the condition p Q p (r) - h(r) ■ 3Q p (r) = -^5p(r)0°(r) 
on the boundary. As Q p = [1 — V]4> p , its perturbative solution shows that it can be viewed 
as a superposition of waves with wavelength ^JD /\ p emanating from a surface localized 
source [1 — V]5p(r). For p — 0, our main focus, the effective source is averaged over a 
diffusion length ^D/Xq, as indicated by the presence of the X p Q p (r) term in its governing 
equation. This leads to Sp 00 ~ as we find in the perturbative solution of the spherical 
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pore|23j and also from exact evaluations |25j. Without any assumptions on the pore geom- 
etry or the p(r) texture, we make some general observations for each contribution to A p . 
Note that the first order term, <5pocb depends sensitively on the symmetry and the profile 
of the mode, (0q) 2 , along the boundary in relation to 5p. While it vanishes for the simple 
situations where 0[] is uniform along the boundary, it may not do so when there exists signif- 
icant variation of 0o as when the complex pore geometry dictates. For empirical Ai (t) with 
a multi-exponential characteristics, as often observed in geophysical applications (26], Eq|3] 
suggests that one cannot safely assume 0o(r) is uniform along the contours of the boundary. 
The chance for a sizeable first order contribution is further enhanced when the texture Sp(r) 
varies commensurate with the former. 

Now we turn to a spherical pore of radius a, and seek exact solutions for both uniform 
Po and p(r) of the form (with a < 1) p(r) = po(l + erf (6)) where we consider the cases of a 
stepwise texture (f(9 < n/2) = —l,f(6 > ir/2) = 1)) and a sinusoidal (f(9) = cos(6 1 )). We 
look for the cigcnmodcs[25j in the form of 0fc(r) = YuL=o s k,LjL(kr)Y®(Q) where is the 

spherical Bessel function, YjP's are the spherical harmonic functions with M = (due to the 
azimuthal symmetry), k represents an infinite set of numbers that allow for a non-trivial 
solution for the coefficients s k that facilitate the boundary condition be met: 

2k^ A LtLl j L (ka)s ktL - (jvika) + ka(j L > +1 (ka) 

L 

-3v-i{ka))s k)L > = (7) 

where £\l,v — J dQf(9)Y®,Y®. Viewed as a homogeneous matrix equation K. ■ s k — 0, the 
eigenvalues are found from the condition that det [K] = 0. We solve this by truncating the 
matrix to a finite though large size and searching numerically for the root. The fractional 
difference in the lowest eigenvalue (Ao = Dk^^) between the uniform and non-uniform cases 
are shown in Figure [2} First panel shows the stepwise texture with a ranging from 0.01 to 
1.0 as indicated. The results from both the exact solution (solid lines) and the second order 
perturbation (points) agree very well for k < 2 for all values of er, while for k > 2, the 
agreement deteriorates progressively as a grows beyond 0.5. Ao itself is shown in the second 
panel. Anomaly occurs in the k ^> 1 limit with o = 1 for which p(r) vanishes on half of 
the hemisphere. In this special case, the limit k ^> 1 acquires a new diffusion controlled 
time scale (i.e. ~ 1/Ao is quadrupled from Aqo = to Aqo/4, as the slowest mode is now 
controlled by diffusion from the p = zone to the other end p — > oo over the distance of 2a. 
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This crossover is missing for the sinusoidal texture with a = 1 (in the third panel) for which 
only a nodal point (8 = n) exists on which p(r) vanishes. This contrasting behavior is also 
verified in numerical simulations. Using the perturbative approach, it is now straightforward 
to incorporate more complicated p(r) textures for small <r, using Figure [2] as a guide for its 
validity. 

Next, we consider the random glass bead pack as an example of realistic porous media for 
which well controlled experiment and simulations could be carried out. Figure [3] shows the 
results from randomwalk simulations |27j based on the Finney pack {28] in which we realize 
three different p textures that clearly violate the conditions necessary for the perturbative 
approach. Type I shows the strongest deviation from the uniform case (and is analogous 
to the case of the hemispherical p of Figj2] with a ~ 1) as we randomly assign a value 
of 0.16 or 1.84 x p to each grain with equal probability. In this case, 5\o/\q/<j 2 ~ 1.22 
is significantly larger than in the closed sphere even though k ~ 0.12. This is likely due 
to the existence of wider spatial separations between the two p values, as expected in the 
pore morphology of a random packing. Type II draws randomly from a distribution of 
p values. Type III uses a texture generated using a correlated random noise[27j. In this 
case, values of Sp(r) are correlated over just a fraction of the bead radius, separating the 
correlations of boundary shape variation from that of Sp(r). Note that even though such 
Sp(r) has a wider distribution (Bottom panels show the histogram and graphical rendition 
of each), the diffusive smearing greatly reduces its effect, and we obtain a result virtually 
indistinguishable from the uniform case. Similar observation had been made numerically 
by Valfouskaya et al[9\. It is beyond the scope of this letter to delve into the relevance of 
these textures for real systems, but they cover a wide range of plausible patterns present in 
natural media such as rocks and certainly realizable in bead packs [29J and other artificial 
structures. [3] Quantitative application of our perturbative solution to realistic 3D systems 
will require the profile of the relevant eigenmodes on the boundary, which may rarely be 
available analytically, but through numerical evaluations of <ftp's. [2"71 130] 
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FIG. 1: Schematics for the difference between the population evolution with uniform po an d an 
inhomogeneous p(r) with the finial slopes given by A[j and Ao as indicated by the broken curves. 
The accompanying change in the spectral weight distribution is reflected in the difference W — W° 
on the y-axis where W = 1 — (ao) 2 and W° = 1 — (a[j) 2 . 
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FIG. 2: Top: Results for a sphere of radius a and varying k = poa/D values with the hemispherical 
p texture. The solid lines show 5Ao/Aq for a = 0.01 ~ 1.0 as obtained from the exact matrix 
formulation described in the text. The filled points represent the second order perturbation result. 

2 

Bottom left panel shows the Ao/Aoo ( Aoo = -Df? ) for varying cr's as above. Bottom-right panel 
compares 5Xq/Xq with a = 1.0 for the hemispherical and the sinusoidal textures of p(r) = po(l + 
a cos (#)). 
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K = 0.12 




FIG. 3: Top: Results for a Finney pack with three differenct p textures (for all poV p /S/D = 0.12) 
as described in the text. Also shown is the result with the matching uniform po for comparison. 
The broken line is an exponential function with the initial slope as expected with X/t s = poS/V p . 
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